#include "PCH.h"
#include "..\Common.h"


AlbersToLatLon::AlbersToLatLon(double p1, double p2, double lt0, double ln0, double S0) {
	// Set up custom region
	phi1=p1;
	phi2=p2;
	lamb0=ln0;
	latorigin=lt0;
	S=S0;

	// Check for a latitude of origin that differs from the midpoint of the two standard parallels
	if (latorigin != (phi1 + phi2)/2) {
		// Find correction factor to adjust for latitude of origin

		dlamb = 0;  // (lamb0 - lamb0)

		// Convert everything to radians for use with sin and cos  
		double dtlat = deg2rad(latorigin);
		double dtdlon=deg2rad(lamb0);
		double tphi1=deg2rad(phi1);
		double tphi2=deg2rad(phi2);
		double tlamb0=deg2rad(lamb0);

		// Perform calculation upon initialization so it's only done once
		rho1 = (2*R*cos(tphi1)) / (sin(tphi1) + sin(tphi2));
		rho2 = (2*R*cos(tphi2)) / (sin(tphi1) + sin(tphi2));

		f1 = (rho1 + rho2) / 2;
		f2 = (sin(tphi1) + sin(tphi2)) / 2;
		f3 = sqrt( (rho1*rho1) + (4*R*R*(sin(tphi1) - sin(dtlat)) / (sin(tphi1) + sin(tphi2)) ) );

		work=dlamb*f2;
		work=deg2rad(work);

		correction = S*(f1-f3*cos(work));
	} else {
		correction = 0;
	}

	// Convert everything to radians for use with sin and cos  
	phi1=deg2rad(phi1);
	phi2=deg2rad(phi2);
	lamb0=deg2rad(lamb0);

	// Perform calculation upon initialization so it's only done once
	rho1 = (2*R*cos(phi1)) / (sin(phi1) + sin(phi2));
	rho2 = (2*R*cos(phi2)) / (sin(phi1) + sin(phi2));

	f1 = (rho1 + rho2) / 2;
	f2 = (sin(phi1) + sin(phi2)) / 2;

}

void AlbersToLatLon::Convert(double x, double y, double *newx, double *newy) {
	double latlon[2] = {0.0};
	// Adjust for the correction factor
	y += correction;
	// Convert to lat/lon
	latlon[0] = (1.0/f2)*atan(x/(S*f1 - y)) + lamb0;
	dlamb= latlon[0]-lamb0;

	double work0=sin(phi1);
	double work1=f2/(2.0*R*R);
	double work2=x/(S*sin(dlamb*f2));

	latlon[1] = asin(work0 - work1*(work2*work2 - rho1*rho1));

	*newy = rad2deg(latlon[0]);
	*newx = rad2deg(latlon[1]);
}